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Abstract:  In  this  study  a  two-dimensional  numerical 
model  for  simulating  ice  transport  and  accumulation 
in  the  vicinity  of  river  ice  booms  is  developed.  The 
model  considers  the  dynamics  of  surface  ice  transport 
in  the  river,  coupled  with  the  hydrodynamics  of 
the  flow.  The  water  flow  inside  the  moving  surface 
ice  and  the  ice  accumulation  is  included  in  the 
hydrodynamics.  The  Lagrangian  discrete-parcel  meth¬ 
od  with  smoothed  particle  hydrodynamics  is  used 
to  simulate  the  ice  dynamics  and  a  finite-element  meth¬ 
od  is  used  to  solve  the  hydrodynamic  equations.  Ice 
entrainment  at  the  boom  or  the  leading  edge  and 


underside  of  the  ice  accumulation,  as  well  as  the  limit¬ 
ing  boom  load  for  ice  retention,  are  considered.  The 
model  is  verified  with  analytical  solutions  for  idealized 
ice  jams  in  a  rectangular  channel,  and  calibrated  to  an 
ice  jam  that  progressed  up  the  lower  Missouri  River 
during  January  1977.  The  model  is  then  used  to  assess 
the  feasibility  of  ice  booms  on  the  lower  Missouri  River. 
The  results  show  that  conventional  ice  booms  may 
not  be  effective  for  typical  flow  conditions  in  the  lower 
Missouri  River,  unless  the  water  level  at  the  Missouri- 
Mississippi  River  confluence  is  high  and  the  water  dis¬ 
charge  is  low. 
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Numerical  Simulation  of 
River  Ice  Control  with  Booms 

LIANU  LIU  AND  HUNG  TAO  SHEN 


INTRODUCTION 

The  middle  Mississippi  River,  from  its  confluence 
with  the  Missouri  River  to  where  it  joins  the  Ohio  River 
at  Cairo,  Illinois  (Fig.  1),  is  a  vital  navigation  route. 
During  the  winter  months,  however,  floating  ice  can 
accumulate  and  block  shipping  from  that  section  of  the 
river.  Ice  sources  are  the  Missouri  River  and  the  middle 
Mississippi  itself.  In  addition  to  suspending  navigation, 
financial  consequences  include  ice  damage  to  river 
training  structures  when  the  ice  releases.  Severe  ice 
conditions  occur  roughly  one  winter  in  seven,  with 
1989,  1979,  1977,  1970,  1962,  1958,  1951,  and  1936 
standing  out  in  the  recent  historical  record.  The  worst 
of  these  cases  occurred  in  January  1977,  when  an  ice 
jam  delayed  barges  at  Cairo  for  27  days  at  an  estimated 
cost  of  $19  million.  The  1977  event  also  caused  $1.5 
million  in  damage  to  river  regulating  structures  below 
Commence,  Missouri,  when  the  ice  jam  released.  Even 
during  winters  without  major  ice  events,  delays  to 
navigation  and  operational  difficulties  resulting  from 
ice  are  common  problems.  Tuthill  and  Mamone  (1998) 
provide  a  detailed  description  of  the  middle  Mississippi 
ice  problems  and  control  alternatives. 

The  Missouri  River,  which  is  officially  closed  in 
winter  for  navigation,  and  is  uncontrolled  for  1250  km 
upstream  from  its  confluence  with  the  Mississippi  River, 
generates  large  amounts  of  ice.  Records  from  the 
Missouri  River  Division  (MRD)  of  the  U.  S .  Army  Corps 
of  Engineers  report  that  ice  jams  have  formed  at  RM* 
6.  The  Lewis  Bridge  at  RM  8  may  aid  the  initial  ice 
arching  (Fig.  2).  The  backwater  from  the  Mississippi 
confluence  may  also  contribute  to  jamming  in  this  reach 
by  reducing  water  surface  slope  and  velocity.  Once 


*River  mile. 


initiated,  this  cover  presents  a  barrier  to  frazil  and  ice 
floes  arriving  from  upstream.  Ice  will  either  accumulate 
at  the  leading  edge  of  the  cover  or,  if  the  water  velocity 
is  sufficiently  high,  it  will  be  carried  under  the  leading 
edge  to  be  deposited  under  the  ice  cover.  If  the  water 
velocity  is  too  high,  ice  will  be  swept  downstream  and 
the  ice  cover  will  not  progress.  Tuthill  and  Mamone 
(1998)  showed  that  installing  an  ice  retention  structure 
in  the  lower  Missouri  River  to  reduce  the  quantity  of 
ice  reaching  the  middle  Mississippi  would  also  have 
the  benefit  of  reducing  the  ice  jam  problems  in  that 
stretch  of  the  river  and  at  the  Missouri-Mississippi 
confluence.  Ice  control  schemes,  consisting  of  floating 
booms,  pier  structures,  and  artificial  islands,  have  been 
used  successfully  on  the  St.  Lawrence  River  (Tuthill 
1995).  For  the  lower  Missouri  River,  the  average  slope 
is  on  the  order  of  1 .7  x  1 0-4.  The  average  water  velocity 
is  high,  with  a  range  of  0.75  to  1.8  m/s.  Using  a  one¬ 
dimensional,  steady-state  model  of  ice  cover 
progression,  Tuthill  and  Mamone  (1998)  found  that,  at 
low  Missouri  River  flows,  the  water  velocity  approaches 
0.7  m/s  between  RM  1 5  and  20,  and  that  floating  booms 
might  be  successful,  so  a  detailed  study  of  the  boom 
locations  was  needed.  This  report  details  a  two- 
dimensional  dynamic  simulation  model  of  river  ice  that 
was  developed  to  assess  the  feasibility  of  ice  control 
on  the  lower  Missouri  River  using  booms. 

MODEL  DEVELOPMENT 

Theories  on  ice  jams  that  are  based  on  the  static 
equilibrium  of  floating  accumulations  of  granular  ice 
are  well  developed  (Pariset  and  Hausser  1961,  Uzuner 
and  Kennedy  1976,  Beltaos  1983).  Static  ice  jam 
theories,  which  neglected  the  dynamic  effect  of  the  ice 


motion,  were  derived  on  the  basis  of  one-dimensional 
formulations.  These  theories  have  been  used 
successfully  to  determine  ice  jam  thickness  along  a  river. 
Flato  and  Gerard  (1986)  and  Beltaos  (1993b)  developed 
one-dimensional  numerical  models  for  computing  the 
configuration  of  static  ice  jams.  Beltaos’  model  is  also 
capable  of  simulating  grounded  jams.  Since  the 
dynamics  of  the  ice  movement  and  flow  were  not 
considered,  the  static  ice  jam  theories  cannot  explain 
the  formation  of  ice  jams.  They  also  cannot  determine 
whether,  when,  and  where  a  jam  will  form.  Moreover, 
the  momentum  effects  of  ice  and  water  flows  on  the  ice 
jam  evolution  and  thickness  were  not  accounted  for  in 
the  static  ice  jam  theories. 

Shen  et  al.  (1990)  developed  an  analytical  frame¬ 
work  for  the  dynamic  transport  of  river  ice  and  ice  jam 
formation  and  evolution.  Lai  and  Shen  (1991) 
developed  a  numerical  model  for  simulating  dynamic 
ice  transport  and  ice  jam  evolution  in  river  channels 
and  showed  the  importance  of  the  inertia  effect  on  ice 
jam  configuration.  They  also  showed  that  the  water 
wave  speed  is  not  affected  by  the  ice  conditions,  and 


that  the  speed  of  the  stress  wave  in  the  ice  layer  is  not 
related  to  the  speed  of  water  waves.  As  a  result  of  this 
independence,  water  and  ice  flow  equations  need  not 
be  solved  simultaneously.  Moreover,  the  speed  of 
characteristic  waves  of  the  water  and  ice  equations  can 
be  significantly  different.  The  temporal  and  spatial 
discretizations  in  the  numerical  solutions  have  to  be 
treated  with  care.  More  recently,  Zufelt  and  Ettema 
(1997)  used  a  simplified  one-dimensional  formulation 
to  study  the  dynamic  effect  on  ice  jam  thickness  profiles 
in  prismatic  channels.  They  used  the  Mohr-Coulomb 
law  for  static  passive  granular  accumulation  to  describe 
the  internal  ice  stress. 

One-dimensional  models  have  limited  applicability, 
as  river  ice  transport  and  ice  jam  evolution  are  two- 
dimensional  phenomena  attributable  to  the  existence 
of  bank  friction  and  nonuniform  water  currents.  Shen 
et  al.  (1993)  developed  a  two-dimensional  dynamic  river 
ice  transport  model,  which  was  successfully  used  to 
study  the  dynamics  of  ice  transport  and  jamming  in  the 
upper  Niagara  River  (Su  et  al.  1997,  Lu  et  al.  1999).  It 
was  also  modified  and  applied  to  study  the  transport  of 
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Figure  2.  Study  area. 


lake  ice  over  the  Lake  Erie-Niagara  River  Ice  Boom 
(Shen  et  al.  1997).  Su  (1997)  further  refined  the  model 
to  include  the  water  flow  in  the  surface  ice  layer,  so 
that  grounded  ice  accumulations  could  be  more 
accurately  modeled. 

In  this  report,  the  models  of  Shen  et  al.  (1997)  and 
Su  (1997)  are  refined  to  simulate  the  dynamic  ice 
transport  and  ice  jam  processes.  This  study  applies  the 
model  to  assess  the  feasibility  of  ice  control  on  the  lower 
Missouri  River  using  floating  booms.  The  model 
equations  and  numerical  methods  are  described  in  the 
following  sections. 


written  as 

aip/r+p/fli-AO+Pifi#]  (i) 

dt 

V(p?i+p7u+pi7icc)  =  0 

in  which 

H'  =  h  +  r ~  water  depth  beneath  the  ice  layer 
h  =  water  depth  below  the  reference  level 
r\ 7  =  elevation  of  the  undersurface  of  ice 
fj  -  ice  layer  thickness 


Hydrodynamic  model 

By  considering  the  flow  in  a  river  with  an  upper 
surface  ice  layer  and  a  lower  water  layer,  as  shown  in 
Figure  3,  the  total  water  and  ice  mass  conservation 
equation  for  the  case  with  a  floating  ice  layer  can  be 


t[  -  Pi^/  =  submerged  ice  thickness 

p  =  water  density 
pi  =  ice  density 

N  =  ice  concentration,  i.e.,  the  volumetric 


3 


Figure  3.  Definition  sketch. 


fraction  of  the  solid  phase  in  the  ice-water 
mixture 

q  j  =  unit-width  water  discharge  beneath  the 
ice  layer 

q  u  —  unit-width  water  discharge  in  the  ice  layer 
q  ice=  Nt-X  -  unit-width  ice  discharge 
V  i  =  ice  velocity. 


Since  the  ice  mass  conservation  gives 


a[£i  tjN] 
dt 


=  -V(Pi?ice) 


eq  1  reduces  to 


#tx  #lx  +  9W  #ty  #ly  #uy  components 
of  total  unit  width  water  discharge 
qlx,  qly  -  components  of  the  unit-width  water 
discharge  beneath  the  ice  layer 
?ux  =  ?ix  +  fe-  tfuy  =  9iy  +  ?sy = water  discharge 
in  the  upper  ice  layer 

<7ix  =  FixOl  —  Tl')(l  -  AO  and  qiy  =  Viy(r\  -110(1 
-  N)  ~  water  discharge  carried  by  ice 
V\x>  v\y  =  ice  velocity 

qsx>  qsy  =  components  of  unit- width  water  discharge 
in  the  ice  layer  relative  to  the  moving  ice 
or  the  seepage  discharge  in  stationary  ice 
accumulations. 


^  +  V-(^+?u)  =  |(^i).  (2) 

Therefore,  the  continuity  equation  for  the  total  water 
discharge  is 


dH  d(qtK)  3(?ty)_  d 

■a7+_a r+~^r~Tt{Nt° 

where 

H  =  h+  ri  =  total  water  depth 
r|  =  water  surface  elevation 


(3) 


When  the  surface  ice  is  grounded,  whether  it  is 
moving  or  stationary,  the  condition  rj  -  rf  =  (pj  /  p)^ 
is  no  longer  valid,  and  the  lower  layer  discharge  q]  is 
zero.  In  this  case,  as  shown  in  Figure  4,  the  water  mass 
conservation  equation  becomes 

Jj[(l  -  A0At|]  =  -V  •  ~q  u  (4) 

in  which,  Ar|  =  H.  Equation  4  can  be  rewritten  as 

dH(l  -  N)  t  dqix  |  3<?ty  0  (5) 

dt  dx  dy 

and  qXx>  qXy  =  0. 


Figure  4.  Definition  sketch  for  grounded  ice  accumulation. 
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The  momentum  equation  in  the  ^-direction  for  the 
under-ice  water  layer  is 

9/  dx  FT  dy  FT  p 


1 

_(^x  + 
p  9jt 


97:x  9n 


and  for  water  in  the  upper  ice  layer  it  is 

dftix  .  j  c^xx,  j  (4ux^uy, 
5/  3;r 

1  Zr^n,  ,/ 

^ix 

n  ox 


Hu  =  (r\-  r|')(l  -  AO  =  net  water  depth  in  the  ice 
layer 

Tix  =  resistance  to  the  upper  layer  flow 
attributable  to  ice 


Tjk  =ejk( 


Mi  ,  a?ik , 


£jk  =  generalized  eddy  viscosity  coefficients 
i  and  k  =  the  two  coordinate  directions 

x  s  =  wind  drag  on  the  water  surface  or  the  ice 
drag  on  the  water  surface  on  the  underside 
of  the  ice  layer 

t  b  -  bed  resistance  to  the  flow 
Mex  -  momentum  exchange  at  the  interface  of  the 
upper  ice  and  lower  water  layer, 

which  is 

_  ,  arr_  K. 

^ex  —  Wwx (  Wwx  0^. 

Wwy  — H-Wwz)  (8) 

where,  vvwx,  and  wwz  are  components  of  water 
velocity  at  the  interface  between  ice  and  water  layers. 
When  eq  6  and  7  are  combined  and  the  unit-width 
discharges  are  expressed  in  terms  of  hydraulic 
conveyance,  i.e.,  q  -  KS p,  the  momentum  equation 
becomes 

dfftx  .  d  ffta  /£l  .^U',  9  gtxgty  | 

9/  dxjcf/r  hJ  a \y  jsf  K/r 

K*  1,  .  (9) 

7f)  =  “(T«+T“-T  bx) 

Hu  p 

dx  p  dx  dy 


K\  -  conveyance  of  the  lower  water  layer 
Ku  =  conveyance  of  the  upper  ice  layer 
Kt  =  total  conveyance 


*U=*8+*i 

Ks  -  conveyance  of  seepage  flow  qs 
Kx  =  conveyance  of  flow  carried  by  ice. 


The  conveyance 

Ks  =  ^  =  Mn-i\') 

where  X  is  a  seepage  coefficient  defined  by  (Bear  1 972) 


X=iki^gds 

in  which 

ds  ” 

p  =  porosity 

Ms  =  average  surface  area  per  volume 
k  =  dimensionless  empirical  coefficient. 

For  randomly  placed  square  plastic  blocks  (5  x  5  x 
0.6  cm  and  10  x  10  x  1.3  cm),  Beltaos  and  Wong  (1986) 
found  k  -  0.70.  Applications  to  the  Credit  River  ice 
jams  by  Beltaos  (1993a)  gave  an  average  X  value  of 
1.6  m/s.  For  the  Restigouche  and  Rushoon  Rivers, 
Beltaos  (1993a)  used  2.5  and  1.0  m/s.  In  this  study  a 
value  of  X  =1.0  m/s  is  used  for  ffeezeup  conditions. 
The  conveyance  K\  is 

K  _  ^(«h  #y_ 

1  % 

The  factor  ah  is  the  fraction  of  the  total  water  flow  depth 
affected  by  the  bed  friction.  Shen  et  al.  (1990)  discussed 
the  distribution  of  shear  stresses  on  the  channel  bed 
and  the  interface  between  the  moving  ice  layer  and  the 
water  current  underneath  it,  and  derived  the  expression 
for  the  coefficient  oeb  as 


nf 


in  which 

Ax  =  flow  area  affected  by  ice  resistance 
Ah  =  flow  area  affected  by  bed  resistance 
=  Manning’s  coefficient  for  ice 
jzb  =  Manning’s  coefficient  for  the  bed 
Fw  =  depth-averaged  current  velocity 
u  =  ice  velocity. 

K-x  is  estimated  as 


|<v*s> 


— ^  ^ 

where  q  j  =  V  i  Hu . 
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Similarly,  the  y-component  of  the  momentum 
equation  is 

dgty  d  9tx<7ty  a  q%  M  Kl 

dt  dx  K2  V  hJ  dy  K\  V  Hj 


-(Tiy+Tsy  ^by) 


(ID 


The  flow  model  solves  for  components  of  q  t  and  the 
water  depth  H  using  eq  3,  9,  and  11.  A  finite-element 
model  with  the  lumping  technique  and  leapfrog  time 
integration  (Connor  and  Brebbia  1978,  Wake  and  Xiao 
1989)  is  used.  The  bed  shear  stresses  can  be  expressed 
as 


?x(<7x+<7  y/2 
Tbx  =  CfP - ~2 - 

<7y(?x+4y)^ 
Tby  =  cfP-^-^ - 


(12) 

(13) 


4-w)=Pcv 

in  which 


Vi-V\ 


(V-Vvy) 


(17) 


j/w  ~-jp=  water  current  velocity 


V  i  =  u  i  +  v  j  =  ice  velocity 

cw  =  water  drag  coefficient  on  ice,  which 
varies  with  ice  concentration  and  ice 
floe  geometry. 

This  coefficient  can  be  related  to  Manning’s  coefficient 
of  the  underside  of  the  cover  as 


°W  -  1/ 

[(1-a  h)/rp 


(18) 


For  a  partially  ice-infested  water  surface,  the  surface 
shear  stress  is  assumed  to  be  a  linear  combination  of 
^a“w)  and  x(‘-w) 


ts  =  (l-A0Ts(a-w,  +  iV  TS 


(i-w) 


(19) 


in  which  the  friction  coefficient  Cf  can  be  expressed  in 
terms  of  Manning’s  coefficients  of  the  bed  and  the  shear 
stress  distribution  coefficient  as 


cf 


ahH" 


g 


On  the  open  water  surface,  the  surface  shear  stress 
attributable  to  wind  effect  can  be  expressed  as 


^rW)=PaY2^2cosea  (14) 

^rW)=PaY2^2sin0a  05) 


in  which 


Y2  =  wind  drag  coefficient  (Wu  1973) 

W  -  wind  velocity  at  10  m  above  the  water 
surface 

pa  =  density  of  air 

0a  =  angle  between  the  wind  direction  and  the 
x-axis. 

For  a  fully  ice-covered  water  surface,  the  surface 
shear  stress  components  on  the  water  can  be  written  as 


The  drag  of  the  seepage  flow  on  ice  is 
isfsx 


"Cix  =  - 


Ki 


^s^sy 


%  =-P^u  2 


(20) 

(21) 


Ice  dynamic  model 

The  momentum  equation  of  the  surface  ice  can  be 
written  in  the  Lagrangian  form  as  (Shen  et  al.  1990) 

— ^ 

F  &+  F  w+  G  (22.) 

Dt 


in  which 
— ^ 

D  V  i  =  acceleration  of  ice 
Dt 

My  =  pjA =  ice  mass  per  unit  area 
R  ~  internal  ice  resistance 


Fa  =  wind  drag 
F  w  =  water  drag 
G  =  gravitational  force 
pj  =  density  of  ice 
N  -  concentration  of  ice 
=  thickness  of  ice. 


4~w)=pcw 


Vi-V , 


(«-*wx) 


(16) 


Force  terms  in  the  momentum  equation  can  be 
expressed  in  two-dimensional  forms  as  follows. 
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(32) 


-> 


Internal  ice  resistance  (R) 

=  Rx  i  +  Ry  j 

*x  =  ^  (axx^i )  “f  (Oxy-^i ) 

(23) 

(24) 

in  which,  <5yy  =  normal  stress  components,  and 

°xy  =  ayx  “  shear  stress  components.  These  stresses 
can  be  determined  by  the  constitutive  relationship. 

— > 

Wind  drag  at  the  air-ice  interface  (F  a) 

A  =  (TaxA0?+(TayA0; 

tax=PaCal^x  (25) 

^=PA\t\Wy  (26) 

in  which 


W  ~J  -  wind  velocity  at  1 0  m  above 

tiie  water  surface 
p  =  density  of  air 
ca  =  wind  drag  coefficient. 

— > 

Water  drag  at  the  ice-water  interface  (F  w) 


—  C^wx^O  1  +  (  ^  wy  N)j 


Uvx  —  Pcv 


fw-fi 


«-^wx) 


(27) 


fw 


-Ki 


(28) 


DA/j 

Dr 


+  MiV-F  = 


0  . 


Since  the  ice  mass  per  unit  area,  Mv  is  determined  by 
the  ice  concentration  and  the  ice  layer  thickness,  one 
more  conservation  equation  is  needed.  The  equation  of 
conservation  of  ice  area  within  an  elemental  area  can 
be  obtained  by  considering  the  ice  area  flux  into  and 
out  of  the  control  area  and  mechanical  redistribution. 


-+NV-V+IL=0  (33) 

Dt 

in  which  Fa  =  rate  of  change  of  ice  area  attributable  to 
mechanical  redistribution. 

A  Lagrangian  discrete-parcel  method  (DPM)  (Shen 
and  Chen  1992,  Shen  et  al.  1993)  is  used  to  simulate 
the  dynamics  of  the  ice  transport.  The  basic  concept  of 
the  discrete-parcel  method  is  that  the  ice,  considered  as 
a  continuum,  can  be  represented  by  a  sufficiently  large 
number  of  individual  parcels.  Each  parcel  has  well- 
defined  properties,  such  as  mass,  concentration, 
thickness,  and  velocity,  and  is  deformable  in  shape.  Ice 
properties  at  parcel  locations  or  finite-element  nodes 
can  be  interpolated  from  the  properties  of  parcels  within 
the  close  vicinity.  The  theoretical  background 
underlying  this  method  is  the  smoothed  particle 
hydrodynamics  (SPH)  developed  by  Lucy  (1977)  and 
Gingold  and  Monaghan  (1977). 

Unlike  the  original  smoothed  particle  hydro¬ 
dynamics,  the  present  discrete-parcel  model  deals  with 
ice  movement  in  a  bounded  domain,  such  as  a  river  or 
a  lake.  A  natural  boundary  condition  at  a  stationary 
boundary  is  a  partial-slip  boundary  condition  with  zero 
normal  ice  flux.  As  an  ice  parcel  moves  along  a  solid 
boundary,  it  is  subjected  to  a  frictional  force.  The 
method  of  images  is  used  in  this  model  for  such  a 
boundary  condition. 

This  model  applies  a  dynamic  Mohr-Coulomb  yield 
criterion  to  calculate  the  boundary  frictional  force  as 
follows 


Gravitational  force  ascribable  to  the  water 
surface  slope  (G) 

(j  =  Gx  ~t +  Gy  ~J 

G^~M'8ltx  (29) 

Gy=~M'8fy-  (30) 

The  ice  mass  conservation  equation  is 

dMj  |  dJI/jU  |  aj/jy_0  (31) 

dt  dx  djy 
This  equation  can  also  be  written  as 


Ff  =  Fz  (34) 

in  which 

Ff  -  frictional  force  between  ice  and  the  solid 
boundary 

Fc  =  ice  cohesive  force,  assumed  to  be  zero 
=  normal  ice  force  against  the  boundary 
<j)B  =  dynamic  friction  angle. 

Constitutive  law 

In  order  to  calculate  the  internal  ice  resistance,  a 
constitutive  law  relating  stresses  with  the  motion  of  ice 
is  required.  The  most  widely  used  constitutive  law  for 
ice  dynamics  is  the  viscous-plastic  law  (Hibler  1979, 
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Wake  and  Rumer  1983,  Shen  et  al.  1993).  The  present 
study  uses  the  Mohr-Coulomb  yield  criterion  (Gutfraind 
and  Savage  1997),  i.e. 

Ci  ±a2  =  (cj  +  a2)sin<|> 

where  Cj  and  c2  are  the  principal  stresses  and  §  is  the 
internal  friction  angle.  The  internal  stress  is  expressed 
in  terms  of  a  nonlinear  viscosity  Jl,  as  follows: 


Cy  =2|iBij-^ekk5ij-P5ij 


(36) 


where  £ij  is  the  strain  rate  tensor.  Assuming  that  the 
principal  axes  of  stress  and  strain  rate  coincide,  we  can 
express  the  viscosity  as 


\i  =  mi  iv 


Psin([) 


>Vn 


ei-£2 


(37) 


along  the  boundary  are  set  to  be  zero.  Since  the 
hydrodynamic  model  is  a  depth-averaged  model,  it  is 
assumed  that  the  floating  ice  boom  does  not  affect  the 
flow  condition.  However,  if  the  boom  configuration  can 
significantly  affect  the  flow,  an  internal  boundary 
condition  may  be  added. 

Ice  dynamics 

Initial  conditions  for  ice  dynamic  simulations  are 
ice  concentration,  ice  velocity,  and  ice  thickness  in  the 
channel.  Ice  thickness,  ice  concentration,  and  ice 
velocity  are  required  upstream  boundary  conditions. 
External  and  internal  forces  acting  on  the  ice  govern 
the  boundary  ice  velocity.  However,  when  the  ice 
concentration  is  not  very  high,  the  internal  ice  resistance 
is  not  significant,  which  is  usually  the  case  at  the 
upstream  boundary,  and  the  ice  velocity  at  the  upstream 
boundary  may  be  approximated  by  the  water  velocity. 


where  ei  and  62  are  the  principal  components  of  the 
strain  rate  tensor  and  fimax,  the  maximum  value  of  the 
viscosity,  determines  whether  the  stress  state  is  inside 
or  on  the  yield  envelope.  When  jl  —  ji^x,  the  ice  will 
flow  as  a  viscous  fluid,  whereas  it  will  flow  in  a  plastic 
manner  when  ji  is  smaller  than  |imax.  The  following 
expression  obtained  by  extending  the  constitutive 
relationship  for  static  ice  jams  (Shen  et  al.  1 990)  is  used 
to  determine  the  pressure  P: 


/>=-(G1+a2)/2  =  ±[l  +  tan2(^±|)] 


L( 


N 


)j 


(38) 


in  which  j  =  1 5,  an  empirical  constant,  and  the  +  and  - 
signs  are  for  convergent  and  divergent  states, 
respectively. 


Initial  and  boundary  conditions 

Hydrodynamics 

Initial  values  of  unit  discharges  qxiqy,  and  water  level 
r|  at  every  finite-element  node  have  to  be  specified  for 
the  hydrodynamic  simulation.  A  steady-state  condition 
that  corresponds  to  the  specified  boundary  conditions 
at  t  =  0  is  generated  with  the  model  and  used  as  the 
initial  conditions  for  qxfqy ,  and  T| .  Boundary  conditions 
used  are  the  water  level  at  the  downstream  boundary 
and  discharge  at  the  upstream  boundary.  The  unit  normal 
discharge  distribution  across  the  upstream  boundary  is 
estimated  by  the  stream-tube  method.  The  calculated 
unit-discharge  distribution  is  adjusted  to  ensure  that  the 
water  level  across  the  width  of  the  boundary  is  constant. 
Along  land  boundaries,  such  as  shorelines  and  dikes, 
the  unit  normal  discharges  at  the  finite-element  nodes 


Limiting  conditions  for  ice  accumulation 
behind  the  boom 

When  surface  ice  arrives  at  the  boom  from  upstream, 
it  may  stop  behind  the  boom  to  accumulate  into  an  ice 
cover.  However,  if  the  current  velocity  exceeds  a  critical 
entrainment  velocity,  the  surface  ice  will  submerge  and 
be  transported  downstream.  If  the  flow  condition 
permits  the  accumulation  of  the  ice  rubble  behind  the 
boom,  the  upstream  progression  is  limited  by  the 
possible  entrainment  of  surface  ice  at  the  upstream  edge 
of  the  ice  accumulation.  In  addition,  an  increase  in  the 
current  velocity  beneath  the  ice  accumulation,  caused 
either  by  the  thickening  of  the  cover  or  an  increase  in 
water  discharge,  can  erode  ice  particles  from  the 
underside  of  the  cover  and  limit  its  progression.  Neither 
surface  ice  entrainment  nor  undercover  erosion  was 
considered  in  the  lake  ice  boom  simulation  model  of 
Shen  et  al.  (1997),  owing  to  the  low  current  velocity  in 
lakes.  The  critical  condition  for  ice  accumulation  behind 
a  lake  ice  boom  is  the  spillover  of  ice  rubble  when  the 
boom  load  exceeds  a  critical  value.  This  condition 
should  also  be  considered  for  river  ice  booms.  However, 
owing  to  the  increased  effect  of  bank  resistance,  the  ice 
load  on  river  ice  booms  rapidly  approaches  a  limiting 
value  as  the  ice  cover  progresses  a  few  channel  widths 
upstream  (Latyshenkov  1946,  Tuthill  and  Gooch  1998). 

In  this  section,  the  limiting  conditions  for  ice 
accumulation  behind  boom  will  be  discussed. 

Ice  entrainment  condition  at  the  boom  and 
the  leading  edge  of  the  ice  cover 

When  water  velocity  is  high,  surface  ice  arriving  at 
the  boom  will  submerge  and  pass  under  it.  Past  field 
experience  suggests  that  ice  retention  is  possible  at  river 
locations  where  the  surface  water  velocity  is  at  or  below 
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0.7  m/s  (2.25  ft/s)  and  the  channel  Froude  number 


does  not  exceed  0.08.  Ashton  (1974)  analyzed  the 
surface  stability  of  floating  ice  blocks  and  obtained  the 
following  submergence  condition: 


Fr  = 


Vc 


2(l4) 


JgtO-y)  t5  H )2^ 


(39) 


where 

Vc  =  water  velocity  upstream  of  the  ice  block 

t b  =  ice  block  thickness 

H  =  flow  depth  upstream  of  the  accumulation. 

This  formula  compared  well  with  experimental  data. 
Larsen  (1975)  studied  the  stability  of  thin  blocks  and 
showed  that  Fc  is  much  larger  than  1.4  when 
(th/H)  <  0.1 .  Daly  and  Axelson  (1990)  further  refined 
these  studies.  More  recent  field  studies  showed  that, 
under  certain  conditions,  specially  designed  booms 
can  perform  successfully  at  water  velocities  as 
high  as  0.76  m/s  (2.5  ft/s)  and  Froude  numbers 
above  0.12  (White  1992).  Results  of  a  recent  1:25 
scale  physical  model  study  by  Tuthill  and  Gooch 
(1998)  support  existing  ice  entrainment  criteria,  find¬ 
ing  that  0.3-m-thick-ice  blocks  submerge  in  a 
4.5-m-deep  channel  at  a  velocity  of  about  0.76  m/s  (2.5 
ft/s)  (all  units  prototype),  and  a  channel  Froude  number 
of  0.1. 


[-2.26(^-)2  + 


2. 14  ^-  +  0.015]^ 


(40) 


where  Vc  =  flow  velocity  below  the  cover  and  L  =  length 
of  ice  floe.  The  experimental  study  of  Kawai  et  al. 
(1997)  showed  that  the  critical  Froude  number  depends 
not  only  on  /L.  For  different  floe  sizes,  their 
experimental  results  showed  that  the  Fcc  =  f(th/L) 
relationships  are  different.  This  shows  the  complexity 
of  the  problem,  and  the  difficulty  in  describing  the 
mechanism  by  a  simple  formula.  For  typical  Missouri 
River  ice  floes  of  thickness  0. 1 5  m,  and  size  L-  1.5 
~  6  m,  eq  40  gives  FCjC=  0.76  ~  1.4  m/s.  In  a  1:25 
physical  model  test  with  natural  ice  scaled  to  the 
observed  ice  size  distribution  of  the  Missouri  River  ice, 
Tuthill  and  Gooch  (1 998)  measured  an  erosion  velocity 
to  be  about  1.5  m/s  (5  ft/s). 

Boom  submergence  condition 

The  present  model  assumes  that  each  span  of  the 
ice  boom  is  designed  to  submerge  and  to  allow  ice  to 
overtop  it  when  a  critical  value  of  cable  tension  is 
exceeded  (Shen  et  al.  1997).  When  the  load  is  reduced 
and  the  weight  of  the  ice  above  the  boom  is  smaller 
than  the  net  buoyancy  of  the  boom,  the  boom  will  rise 
again  to  prevent  ice  from  passing  through.  The  value 
of  critical  cable  tension,  which  varies  with  the  geometry 
and  size  of  booms,  has  to  be  specified  for  the  model 
simulation. 


Erosion  of  ice  on  the  underside  of  ice  jams 

When  an  ice  jam  forms,  its  thickness  may  be  limited 
by  the  stability  of  ice  particles  on  its  underside.  As  the 
jam  thickens  near  its  downstream  end,  the  water  flow 
area  decreases,  increasing  water  velocity  and  shear 
on  the  ice  underside.  The  shear  may  become  large 
enough  to  erode  ice  pieces  from  the  jam’s  underside 
and  transport  them  downstream.  Thinning  the 
downstream  end  of  the  jam  by  erosion  lowers  the  water 
level  at  the  upstream  end,  increasing  the  velocity  and 
the  tendency  for  ice  entrainment.  This  process  of  under¬ 
ice  erosion  near  the  toe  and  entrainment  at  the  upstream 
end  of  the  jam  may  be  sufficient  to  halt  upstream 
progression. 

The  stability  of  ice  floes  on  the  underside  of  an  ice 
cover  or  ice  jam  has  been  investigated  by  Ashton  (1974), 
Uzuner  (1977),  Tatinclaux  and  Gogus  (1981),  and  Daly 
and  Axelson  (1990).  Tatinclaux  and  Gogus  (1981) 
recommended  an  empirical  stability  criterion  in  terms 
of  the  critical  Froude  number  Fe  c  and  the  floe  aspect 
ratio  rb  /L 


Ice  load  on  the  boom 

When  the  boom  stops  the  ice  movement,  mechanical 
thickening  will  occur  and  the  internal  stress  will  build 
up  quickly  in  the  ice  rubble.  The  load  on  the  boom  and 
the  span  cable  tension  can  be  calculated  from  the  ice 
rubble  stresses. 

Consider  a  triangular  differential  element  of  ice 
rubble  abc  (Fig.  5)  in  contact  with  the  boom,  which  has 
a  finite  thickness  tv  The  total  force  acting  in  the  x- 
direction  is 

£FX  =  -GxxNt\dy  -  Gy^Nt^dx  +  XNt\dl  + 
j  lsxNdxdy  +  ^  zwxNdxdy  +  j  p  igNt-ldxdy ^  ' 

(41) 

The  total  force  acting  in  the  ^-direction  is 

£Fy  =  -OyyNt^dx  -  Oy-yNt^dy  +  YNt^dl  + 
j  xsyNdxdy  +  j  xwyNdxdy  +  ±p{  gN^dxdy^ 

(42) 
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X 

Figure  5.  Interaction  of  a  differential  element  of  ice  rubble  with  the  boom. 


in  which 

X  and  Y  =  x-  and  ^-components  of  the  ice  load  on  the 
boom 

N  =  ice  concentration 

Tsx,Tsy  =  components  of  wind  drag  in  x-  and 
^-directions 

Twx’Twy  =  components  of  water  drag  in  x-  and 
y-directions. 

Since  the  ice  element  is  essentially  stationary,  the  force 
balance  conditions  X  =  0  and  X  “  0  reduce  the 
above  equations  to  the  following,  by  neglecting  the 
higher  order  terms 


axx  cos(//,.r)  +  ayx  cos (n,y) 

(43) 

Y=  Oyy  cosO,  y) + cxy  cosO,  x)  . 

(44) 

The  normal  and  shear  stresses  acting  on  the  boom 
segment  ab  are 

Gn  =  oxx  cos20,  x)  +  Gyy  cos20,  y)  + 

(45) 

2axy  cos(w,  x)  cos  (n,  y) 

ct  =(oyy  ~Gxx)cos(«,x)cos(rt,y)  + 

(46) 
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Gxycos  (n,x)-ayx  cos  (n,y). 

The  normal  and  tangential  components  of  the  load  per 
unit  length  of  the  boom  are 

Pn=Mian 

(47) 

Px=NtxGi  . 

(48) 

MODEL  VERIFICATION  AND  CALIBRATION 

Verification  of  the  hydrodynamic  model  with 
a  uniform  ice  cover 

The  hydrodynamics  model  is  verified  by  simulating 
flow  in  a  rectangular  ice-covered  channel  of  length  2000 
m,  width  500  m,  bed  elevation  at  upstream  5  m,  bottom 
slope  0.0001,  and  bed  Manning’s  coefficient  of  0.03. 
The  ice  cover  thickness  is  selected  to  be  1 .06  m,  with  a 
Manning’s  coefficient  0.03.  Water  discharge  is  2500 
m3/s  and  downstream  water  surface  elevation  is  -0.2 
m.  The  simulated  water  surface  profile  is  compared  with 
the  result  of  the  one-dimensional  backwater  calculation 
in  Figure  6. 

Analytical  solutions  for  river  ice  jams 

Analytical  solutions  for  idealized  ice  jams  in  a 
straight,  uniform,  rectangular  channel  with  a  uniform 
current  are  developed  for  verifying  the  numerical  model. 
The  numerical  simulation  results  are  compared  with 
analytical  solutions  for  ice  accumulation  in  a  channel 
of  length  5000  m,  width  500  m,  constant  water  velocity 
of  0.6  m/s,  and  zero  water  surface  slope.  A  boom  is 
placed  at  500  m  and  assumed  100%  effective,  i.e.,  no 
ice  is  allowed  to  pass  it.  Initially,  900  ice  parcels  of  size 
50  x  50  m  having  a  thickness  of  0.2  m  and  a 
concentration  of  0.6  are  placed  over  the  water  surface 
from  the  upstream  boundary  and  the  boom.  In  the 
simulations,  this  ice  cover  was  allowed  to  move  and 
consolidate  behind  the  boom  under  the  action  of  the 
current  drag.  No  additional  ice  was  added  during  the 
simulations.  The  model  parameters  used  are 
summarized  in  Table  1 . 
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Figure  6.  Comparison  of  water  surface  profiles. 


Without  bank  friction 

The  steady-state  ice  thickness  profile  can  be  obtained 
by  simplifying  the  momentum  equation,  eq  22.  When 
only  the  water  drag  inx-direction,  i.e.,/^  =  P^w^wx  > 
is  considered  and  the  bank  friction  is  neglected,  the  ice 
momentum  equation  is  simplified  to  Rx  +  F ^  =  0.  The 
internal  ice  resistance  reduces  to 

A  =  £(°xx^i ) +  )  =  Yx(PMi ]  ■ (49) 


With  bank  friction 

An  analytical  solution  for  the  width-averaged  ice 
jam  thickness  profile  can  be  derived  by  extending  the 
solution  of  Pariset  and  Hausser  (1961). 

2g]  J_ 

/i=4qO-*-M2  (51> 

where  B  is  the  channel  width  and  the  equilibrium  ice 
thickness  is 


Using  eq  37  for  the  pressure  term  leads  to  a  simple 
analytical  solution  for  the  static  ice  accumulation 
thickness  profile 

h = (4 + — ^5^ —  (so) 

tan2(f  +  -f)(l-^)Piir 

in  which  £i0  =  single  layer  ice  thickness  and  x  -  distance 
from  the  leading  edge  of  the  jam  where  q  =  t-^.  The 
simulated  results  and  analytical  solution  are  compared 
in  Figure  7.  Figure  7a  shows  progressive  thickening 
and  compressing  of  the  simulated  ice  cover  with  time. 
Figure  7b  compares  the  analytical  solution  to  the 
simulated  ice  thickness  profile  on  and  after  t  =  4  hours. 


4q  -  (- 


BNCJ\ 


■)2 


where 

|i2  =  N  tan<()(l  +  sin4>) 

and 

Pj  =  tan(})(l  —  sin(J))  (Beltaos  1995) 
For 


(J>  =  46° 
p2  =  1*068 
Pj  -  0.29. 


(52) 


Table  1. 

Parameters  used  in  the  ice  dynamic 

simulation. 

Parameter  Description 

Value 

Maximum  ice  concentration 

0.6 

* 

Internal  friction  angle  of  ice 

46° 

tan  (|) 

Boundary  friction  coefficient 

1.04 

j 

Empirical  constant 

15 

cw 

Water  drag  coefficient  on  ice 

0.02 
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a.  Simulated  time-dependent  evolution  of  ice  jam  profile. 


Distance  Along  River  (m) 

b.  Comparison  of  final  jam  profile  with  steady-state  analytical  solution. 

Figure  7.  Comparison  of  simulated  jam  profiles  and  the  analytical  solution 
for  the  zero  bank  friction  case. 

The  simulated  result  and  analytical  solution  are  In  this  study,  the  constitutive  law  is  modified  for  a  small 
compared  in  Figure  8.  The  simulated  profile  for  t  >  4  strain  rate  condition  to  avoid  these  problems, 

hours  and  the  analytical  solution  are  different  because  In  the  simulation,  critical  values  of  and  P  for 

the  ice  momentum  was  neglected  in  the  analytical  each  ice  parcel  were  determined  as  the  following 

solution.  Figure  9  presents  two-dimensional  plots  conditions  are  all  satisfied:  1)  ice  parcel  velocity  is 

showing  the  development  of  the  ice  thickness  smaller  than  a  small  critical  value  of  0.001  m/s,  2)  ice 

distribution  upstream  of  the  boom.  parcel  velocity  is  smaller  than  the  value  in  the  previous 

Theoretically,  the  viscous-plastic  constitutive  law  time  step,  and  3)  5=|ei-B2  I  is  smaller  than  a  very 

cannot  simulate  static  conditions  correctly.  For  very  small  value  5C  =  1  x  10"^s_f  When  these  critical 

small  strain  rates,  the  viscosity  becomes  very  large.  In  conditions  are  reached,  the  following  approximation  is 

numerical  computations,  a  limiting  value  for  \l  is  often  introduced  to  calculate  ice  stresses  for  the  particular 

used.  This  approximation  changes  the  constitutive  ice  parcel: 
relationship  to  a  linear  viscous  law  and  the  shear  stress  p 

approaches  zero  when  the  strain  rate  approaches  zero.  Gij =  aijc  “Tf  (53) 
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a.  Simulated  time-dependent  evolution  of  ice  jam  profile. 


Figure  8.  Comparison  of  simulated  jam  profiles  with  the  analytical  solution 
for  the  case  with  bank  friction. 


where  ayc  and  Pc  are  critical  values  of  ay  and  P, 
respectively.  In  addition,  when  the  velocity  of  an  ice 
parcel  is  less  than  0.5  mm/s,  the  parcel  is  stopped.  Using 
this  method,  we  can  simulated  the  static  ice  jam 
condition  with  the  viscous-plastic  constitutive  law. 

Open  water  calibration  for  Missouri  River 

The  model  is  calibrated  for  the  open  water  condition 
to  determine  the  Manning’s  coefficient  of  the  channel 
bed  so  that  the  model  can  correctly  simulate  the  current 
velocity  distribution  and  water  surface  slope.  The 
detailed  geometry  survey  data  provided  by  Kenneth 
Balk  &  Associates,  Inc.,  cover  the  Missouri  River  from 


RM  0  to  24,  as  well  as  the  Mississippi  River  from  St. 
Louis  to  the  mouth  of  the  Missouri  (RM  1 87-195),  with 
cross  sections  about  300  m  apart.  Tuthill*  developed 
the  stage-discharge  relationships  using  HEC-2 
simulations  for  the  reach  from  the  St.  Louis  gage  up  to 
about  RM  30  on  the  Missouri  River  (see  Fig.  1  for  a 
map  of  the  confluence  area).  The  HEC-2  model  was 
calibrated  to  flows  in  September  1994  at  Hennann, 
Missouri,  the  corresponding  observed  water  surface 
elevations,  and  average  velocity  from  USGS  stream 


*Personal  communication  with  A.  Tuthill,  CRREL,  1997. 
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b.  Hour  4. 

Figure  9.  Ice  thickness  distributions. 


gaging  data.  The  present  hydrodynamic  model  is 
calibrated  to  the  HEC-2  results. 

In  a  one-dimensional  study  of  ice  jamming  on  the 
Missouri  using  the  ICETHK  model,  Tuthill  found  two 
possible  locations  where  an  ice  boom  may  be  used:  RM 
1 6  and  RM  8.2.  For  computational  efficiency,  two  model 
domains  are  used  for  these  two  boom  sites.  Domain  1 
covers  RM  13-20  and  domain  2  covers  RM  5-13.  The 
finite-element  mesh  and  bed  elevation  for  domain  1  are 
shown  in  Figure  10.  Similar  plots  for  domain  2  are 
shown  in  Figure  1 1 .  The  two-dimensional  model  water 
surface  profiles  are  compared  to  HEC-2  water  levels  in 
Figure  12.  The  Manning’s  coefficients  for  the  bed  are 
shown  in  Table  2. 

Simulation  of  the  January  1977  Missouri  River 
Ice  Jam 

Manning’s  coefficients  for  the  river  bed  were 
determined  by  the  open  water  calibration  in  the  last 
section.  The  Manning’s  coefficient  for  the  ice  cover  had 
to  be  calibrated,  and  the  internal  friction  angle  of  ice  <j) 


and  the  empirical  constant  j  in  eq  37  had  to  be  validated. 
The  MRD  ice  data  spread  sheets  for  the  January  1977 
ice  event  give  the  daily  position  of  the  leading  edge  of 
a  ffeezeup  ice  jam  that  progressed  upstream  during  1 7- 
22  January  1977.  Using  an  average  flow  of  650  m3/s 
(23,000  ft3/s),  (j)  =  45°,  jn2  =  1.17,  and  =  0.04,  the 
ICETHK  model  simulation,  with  the  ice  jam  toe  placed 
below  the  confluence  at  RM  193.8  of  the  Mississippi 
River,  gives  an  average  ice  jam  thickness  of  about  0.9 
m  (3  ft)  upstream  of  RM  8  on  the  Missouri  River.*  The 
ICETHK  model  is  a  one-dimensional  steady-state 
model,  which  calculates  the  ice  thickness  profile  of  a 
wide-river  jam  based  on  an  ice  jam  force  balance  that 
is  similar  to  the  model  of  Flato  and  Gerard  (1986). 
Assuming  an  ice  jam  porosity  of  0.4,  we  determined 
the  average  ice  discharge  upstream  of  the  jam  to  be 
about  13.5  m3/s  (475  ft3/s)  to  account  for  this  observed 
upstream  progression.  Assuming  a  floe  thickness  of  0. 1 5 


*Personal  communication  with  A.  Tuthill,  CRREL,  1997. 
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a.  Finite-element  mesh  for  RM  13-20. 


b.  Bed  elevation  map  for  RM  13-20. 
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c.  Finite-element  mesh  for  RM  16-18. 
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b.  Bed  elevation  map  RM  5-13. 
Figure  11.  Model  domain  2. 
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d.  Bed  elevation  for  RM  8-10. 
Figure  11  (conf'dj. 
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a.  Domain  1  (RM  13-20). 


b.  Domain  2  (RM  5-13). 

Figure  12.  Comparison  of  ice-free  water  surface  profiles. 


m  (0.5  ft)  and  an  ice  velocity  of  1.2  m/s  (4  ft/s),  we 
estimated  the  surface  concentration  of  the  moving  pans 
and  floes  to  be  about  30%.  This  estimate  compares 
reasonably  well  with  aerial  photographs  of  moving  ice 
at  the  mouth  of  the  Missouri  on  9  January  1979,  9 


Table  2.  Calibrated  bed  roughness. 

RM 

Manning’s  n  for  bed 

5-7 

0.022 

7-9 

0.020 

9-11 

0.015 

11-13 

0.015 

13-20 

0.019 

February  1989,  and  16  December  1989.  These  photos 
showed  surface  ice  concentrations  of  20  to  30%.*  On 
the  basis  of  this  ice  discharge  estimate  and  a  downstream 
water  surface  elevation  estimated  by  the  ICETHK 
model,  the  ice  accumulations  behind  the  booms  at  RM 
16  in  domain  1  and  RM  8.2  in  domain  2  were  simulated 
by  the  present  two-dimensional  model. 

Simulation  with  a  boom  at  RM  16 

The  downstream  boundary  water  level  for  domain 
1  at  RM  1 3  was  obtained  by  making  the  simulated  water 
level  at  boom  location  RM  1 6  the  same  as  the  ICETHK 
model  result,  which  was  128  m  (420  ft).  The  Manning’s 


♦Personal  communication  with  A.  Tuthill,  CRREL,  1997. 
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c.  Hour  20. 

Figure  13.  Simulated  ice  thickness  distribution  behind  the  boom  at  RM  16. 


coefficient  n \  is  0.02  for  single  ice  layer  and  increases 
with  the  ice  thickness  (Shen  and  Chen  1992).  The  nx 
value  was  limited  to  a  maximum  of  0.05  for  multilayer 
ice.  Water  discharge  was  650  m3/s  (23,000  ft3/s).  At 
the  upstream  boundary,  the  ice  concentration  is  set  to 
be  0.45  in  the  main  stream  of  the  river,  which  gave  an 
ice  discharge  of  about  7.4  m3/s  (260  ft3/s)  in  the 
simulation. 

The  simulated  ice  thickness  distributions  at  hours 
10,  15,  and  20  are  shown  in  Figure  13.  At  hour  20  the 
simulated  ice  cover  was  about  3.2  km  (2  miles)  long, 
with  a  thickness  in  the  0.6-  to  0.9-m  (2-  to  3-ft)  range. 


In  this  case,  the  water  level  at  the  boom  location 
was  very  high  because  of  the  assumption  that  an  ice 
cover  existed  downstream  of  RM  13  at  the  start  of  the 
simulation.  The  simulated  maximum  water  velocity  near 
the  boom  was  only  0.6  m/s  (2  ft/s),  and  no  ice  passed 
the  structure. 

Simulation  with  a  boom  at  RM  8.2 

In  the  simulation  done  for  domain  2  (RM  5-13), 
with  a  boom  at  RM  8 .2,  the  downstream  boundary  water 
level  at  RM  5  of  124.36  m  (408  ft)  was  obtained  from 
the  ICETHK  simulated  ice  jam  profile  initiated  at  the 


a.  Hour  10. 


b.  Hour  15. 

Figure  14.  Simulated  ice  thickness  distribution  behind  the  boom  at  RM  8.2. 
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c.  Hour  20. 
Figure  14  (cont’d). 


confluence.  The  water  level  at  the  boom  location  is 
125.27  m  (41 1  ft).  Water  discharge  is  650  m3/s  (23,000 
ft3/s),  the  same  as  in  the  previous  case.  Ice  discharge  is 
about  6.8  m3/s  (240  ft3/s). 

Since  the  water  velocity  near  the  boom  is  less  than 
0.6  m/s  (2  ft/s),  ice  stops  at  the  boom,  and  an  ice 
accumulation  develops  upstream. 

The  simulated  ice  thickness  distributions  at  hours 
1 0, 1 5,  and  20  are  shown  in  Figure  14.  The  results  show 
that  the  ice  cover  length  is  about  3.2  km  (2  miles)  at 
hour  20,  and  the  thickness  was  in  the  range  of  0.9  to 
1 .2  m  (3  to  4  ft),  similar  to  the  ICETHK-calculated  ice 
thickness. 


ASSESSMENT  OF  ICE  BOOM 
ALTERNATIVES 

The  calibrated  model  was  used  to  further  assess  the 
effectiveness  of  ice  boom  alternatives.  Simulation 
results  are  presented  here. 

100%  Effective  ice  booms  normal  to 
flow  direction 

Two  simulations  were  carried  out  with  a  100% 
effective  ice  boom.  Limiting  conditions  for  ice 
accumulation  behind  the  boom  were  not  imposed  in 
these  simulations.  The  initial  condition  in  the  first 
simulation  for  domain  1  (RM  13-20),  with  the  boom  at 
RM  16,  was  an  open  water,  steady-state,  ice-free  flow 
at  a  discharge  of  566  m3/s  (20,000  ft3/s).  The 


downstream  boundary  condition  at  RM  1 3  was  a  water 
surface  elevation  of  123.6  m  (405.6  ft).  Upstream 
boundary  conditions  were  a  constant  water  discharge 
of  566  m3/s  (20,000  ft3/s)  and  an  ice  discharge  of  1 1.3 
m3/s  (400  ft3/s).  In  the  second  simulation,  for  domain  2 
(RM  5-13),  the  boom  was  placed  at  RM  8.2,  above  the 
Lewis  Bridge.  The  same  water  and  ice  discharge 
conditions  as  in  the  first  simulation  were  used.  The 
downstream  water  level  at  RM  5  was  120.95  m  (396.8 
ft).  The  goals  of  these  simulations  were  to  see  how  the 
ice  jams  and  boom  loads  would  develop  if  the  booms 
were  100%  effective. 

The  simulated  results  showed  that  the  ice  thickness 
at  the  boom  locations  was  high  and  that  the  ice 
grounded  across  most  of  the  channel  width  near  the 
boom,  as  shown  in  Figures  15  and  16.  The  simulated 
jams  increased  upstream  water  levels  significantly,  as 
shown  in  Figure  17.  Upstream  of  the  toe  area,  the  ice 
thickness  decreased  to  about  0.9  to  1 .2  m  for  the  boom 
located  at  RM  8.2,  and  0.9  to  1 .8  m  for  the  boom  located 
at  RM  1 6. 

The  loads  on  the  booms  are  shown  in  Figure  18. 
The  maximum  load  per  unit  width  was  about  40  to  50 
kN/m.  The  results  showed  that  the  boom  loads  leveled 
off  as  the  jam  extended  upstream.  This  was  attributable 
to  the  grounding  of  the  jam,  as  well  as  the  increase  in 
bank  resistance  as  the  jam  progressed  further  upstream. 
In  domain  1,  because  of  the  low  water  discharge,  the 
high  bed  elevation  causes  the  right  one-third  of  the 
channel  width  to  have  no  flow. 
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Figure  15.  Simulated  ice  thickness  distribution  behind  the  100%  effective  boom  at  RM  16. 


Figure  16.  Simulated  ice  thickness  distribution  behind  the  100%  effective  boom  at  RM  8.2. 
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a.  For  the  boom  at  RM  16. 


b.  For  the  boom  at  RM  8.2. 
Figure  17.  Water  surface  profiles. 


Slant-oriented  boom 

As  discussed  earlier,  when  the  open  water  velocity 
approaches  0.7  m/s,  the  ice  pans  and  floes  approaching 
the  boom  will  underturn  and  flow  downstream, 
rendering  the  boom  ineffective.  Upstream  of  the 
backwater  effect  of  the  confluence,  the  water  velocity 
in  the  lower  Missouri  River  is  typically  above  the 
threshold  for  ice  retention  by  conventional  means.  Even 
at  the  extremely  low  flow  of  5 66  m3/s  (20,000  ft3/s), 


velocity  is  in  the  range  of  0.75-0.9  m/s.  Specially 
designed  booms  may  be  able  to  form  an  ice  cover  under 
these  high  velocity  conditions,  however.  For  example, 
a  simulation  with  a  slant-oriented  boom,  instead  of  a 
boom  normal  to  the  flow,  is  shown  in  Figure  19.  The 
boom  is  placed  at  RM  1 6  at  an  angle  of  about  45°  to  the 
flow.  This  alignment  reduces  the  velocity  component 
normal  to  the  boom  to  70%  of  the  longitudinal  velocity 
component. 
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a.  At  RM  16. 


b.  At  RM  8.2. 


Figure  18.  Simulated  load  distribution  on  the  boom. 


In  the  simulation,  as  the  ice  parcels  reached  the 
boom,  they  were  not  entrained  underneath  but  moved 
along  the  boom  to  accumulate  in  the  shallow  area  near 
the  shore,  initiating  progression  upstream.  With  the 
increase  in  ice  jam  thickness,  the  water  velocity 
underneath  the  ice  cover  increases  as  well.  The  critical 
erosion  velocity  of  1 .5  m/s  was  exceeded  at  the  toe  of 
the  jam  and  ice  started  to  be  transported  downstream. 
This  prevented  the  cover  from  progressing  further 
upstream. 


SUMMARY 

In  this  study  a  two-dimensional  numerical  model 
for  simulating  ice  transport  and  accumulation  behind 
river  ice  booms  was  developed.  The  model  considered 
the  dynamics  of  surface  ice  transport  in  the  river, 
coupled  with  the  hydrodynamics  of  the  flow.  The 
Lagrangian  discrete-parcel  method  with  smoothed 
particle  hydrodynamics  was  used  to  simulate  the  ice 
dynamics,  and  a  finite-element  method  was  used  to 


25 


Bed  Elevation  (m) 


300  m 

i _ i 


Figure  19.  Bed  elevation  in  the  vicinity  of  the  slant-oriented  boom  at  RM  16. 


solve  the  hydrodynamic  equations.  Ice  entrainment  at 
the  boom  and  the  leading  edge  of  the  ice  accumulation, 
erosion  of  ice  on  the  underside  of  the  ice  accumulation, 
and  the  limiting  ice  boom  load  for  ice  retention  were 
considered.  The  model  was  verified  with  analytical 
solutions  for  idealized  ice  jams  in  rectangular  channels, 
and  calibrated  to  an  ice  jam  that  progressed  up  into  the 
lower  Missouri  River  from  the  middle  Mississippi  River 
during  January  1977.  The  calibrated  model  was  then 


used  to  study  the  feasibility  of  using  ice  booms  to  retain 
the  ice  in  the  lower  Missouri  River  to  reduce  the 
jamming  potential  in  the  Mississippi  River  at  the 
Missouri-Mississippi  confluence  and  the  middle 
Mississippi  River. 

Numerical  simulations  were  made  with  a  boom 
located  at  RM  16  or  RM  8.2,  the  two  most  favorable 
locations  for  an  ice  boom.  Table  3  summarizes  the 
simulation  results.  In  the  first  group  of  simulations, 


Table  3.  Summary  of  simulation  results. 

Simulation 

cases 

Water 

discharge 

(m3/s) 

Ice 

discharge 

(m3/s) 

Model 

domain 

(RM) 

Water  level  at 
downstream 
boundary 

Boom 

location 

(RM) 

Comments 

1.  Simulation  for 

Jan.  1977  jam 
event 

650 

7.36 

13-20 

127.10 

16 

High  downstream  water  surface 
elevation  attributable  to  the 
presence  of  downstream  ice  jam. 

650 

6.79 

5-13 

124.36 

8.2 

Boom  is  effective,  fj  =  0.6-1 .2  m, 
similar  to  ICETHK  result. 

2.  100%  effective 
booms 

566 

11.32 

13-20 

123.63 

16 

With  normal  open  water 
conditions,  downstream  water 
levels  are  much  lower  than  in  the 
two  cases  above. 

566 

11.32 

5-13 

120.94 

8.2 

Ice  partially  grounded  near  the 
boom.  Boom  loads  as  high  as  40 
to  50  kN/m. 

3.  Boom  oriented 
at  45°  to  the  flow 

566 

11.32 

13-20 

123.63 

16 

The  boom  may  stop  ice  floes,  but 
erosion  and  entrainment  limit 
upstream  progression. 
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water  levels  at  the  downstream  boundary  were  assumed 
to  be  the  same  as  those  during  the  January  1977  ice 
jam,  which  extended  from  RM  193.8  of  the  Mississippi 
River  to  upstream  of  the  boom  locations.  The  ice  and 
water  discharges  of  the  January  1977  ice  jam  were  also 
used.  These  simulations  showed  that  the  boom  could 
be  effective  at  both  locations.  In  the  second  group  of 
simulations,  an  initial  open  water  condition,  with  a  low 
river  discharge  of  566  m3/s,  was  assumed.  The  ice 
booms  were  assumed  to  be  100%  effective,  i.e.,  the 
limiting  conditions  for  ice  accumulation  behind  the 
boom  were  not  imposed.  These  simulations  showed  that 
the  lower  downstream  water  levels  of  the  open  water 
condition  resulted  in  partial  grounding  of  ice  in  the 
vicinity  of  the  booms.  The  boom  loads  were  not  uniform 
across  the  width,  and  leveled  off  as  the  ice  cover 
extended  upstream.  The  ice  loads  on  the  booms  reached 
very  high  values,  ranging  from  40  to  50  kN/m.  Since 
the  results  of  the  second  group  of  simulations  showed 
that  the  flow  velocity  could  exceed  the  entrainment 
velocity,  an  additional  simulation  was  made  with  a  boom 
oriented  at  an  angle  of  45°  to  the  main  flow  direction. 
This  reduced  the  normal  velocity  of  ice  floes  approaching 
the  boom,  and  the  ice  was  stopped  at  the  boom.  However, 
as  the  ice  accumulation  thickened  behind  the  boom,  the 
critical  undercover  erosion  velocity  was  exceeded  and 
the  cover  progression  was  halted. 

Before  the  present  study,  ice  retention  behind  booms 
was  considered  to  be  marginally  possible,  based  on 
existing  critical  water  velocity  and  Froude  number 
criteria.  The  one-dimensional  steady-state  ice  jam  model 
ICETHK  predicted  that  a  stable  ice  accumulation  was 
possible,  but  the  model  was  unable  to  address  dynamic 
processes.  During  problem  ice  years  on  the  middle 
Mississippi  River,  the  average  Missouri  River  discharge 
is  about  850  m3/s  (30,000  ft3/s).  Two-dimensional 
simulations  at  the  two  most  favorable  sites  (RM  8.2 
and  1 6)  carried  out  here  found  that  ice  retention  behind 
booms  is  unfeasible  even  at  a  much  lower  discharge  of 
566  m3/s  (20,000  ft3/s),  unless  a  specially  designed 
boom  with  high  ice  retention  capacity  can  be  developed. 
The  two-dimensional  dynamic  ice  transport  model 
proved  to  be  a  valuable  tool  for  addressing  important 
design  issues  that  could  not  be  answered  by 
conventional  methods. 
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